%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time-dependent optical Bloch equations for a 3-level system with 2 lasers
%  in V-configuration
%  n = [11 22 33 12 21 13 31 23 32]
%
%  --- |1>
%    \  --- |3>	
%	  \ /
%     --- |2>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Bloch

    o12=	      1;                                                        % Rabi frequencies
    O23=[0:.02:1]*5*o12;
    G12=	.01*o12;                                                        % decay rates
    G23=	    G12;
    G13=   1e-4*G12;
    d12=          0;
    d23=          0;

    n0 =[0 1 0 0 0 0 0 0 0]';
    Tim=[0:.005:1]*10;
    for jj=1:length(O23)
        for tt=1:length(Tim)
            t  =Tim(tt);
            o23=O23(jj);                                                    % detunings
            n  = Matrix_V(n0,t,G12,G23,G13,o12,o23,d12,d23);
            nt(jj,tt)=real(n(1));
            nc(jj,tt)=imag(n(4));
        end
    end
    
    figure(1);      clf;
    subplot(221),mesh(Tim,O23,nt)
    set(gca,'FontSize',16)
    xlabel('\fontsize{16}{\itt}')
    ylabel('\fontsize{16}\Omega_{23}')
    zlabel('\fontsize{16}{\it\rho}_{11}')
    text(12,5,1,'\fontsize{16}(a)')
    shading flat
    axis tight
    view([190 50])
    subplot(222),mesh(Tim,O23,nc)
    set(gca,'FontSize',16)
    xlabel('\fontsize{16}{\itt}')
    ylabel('\fontsize{16}\Omega_{23}')
    zlabel('\fontsize{16}Im {\it\rho}_{12}')
    text(12,5,1,'\fontsize{16}(b)')
    shading flat
    axis tight
    view([190 50])
end


function n = Matrix_V(n0,t,G12,G23,G13,o12,o23,d12,d23)
    I=          1i/2;
    g12=     G12+G13;
    g13= G12+G23+G13;                                                       % decay of coherences
    g23=         G23;
    d13=     d12-d23;
    L12=g12/2+1i*d12;
    L13=g13/2+1i*d13;
    L23=g23/2+1i*d23;
    
% M=[	0	g12	g13	I*o12	-I*o12	I*o13	-I*o13	0	0
% 	0      -g12-g23 0	-I*o12	I*o12	0	0	0	0
% 	0	g23	-g13	0	0	-I*o13	I*o13	0	0
% 
% 	I*o12	-I*o12	0	-L12'	0	0	0	I*o13	0
% 	-I*o12	I*o12	0	0	-L12	0	0	0	I*o13
% 
% 	I*o13	0	-I*o13	0	0	-L13'	0	I*o12	0
% 	-I*o13	0	I*o13	0	0	0	-L13	0	I*o12
% 
% 	0	0	0	I*o13	0	I*o12	0	-L23'	0
% 	0	0	0	0	I*o13	0	I*o12	0	0  ];
    M=[
       -G12-G23 0        G13	 I*o12	-I*o12	0		0		0		0
        G12     0        0		-I*o12	 I*o12	0		0		I*o23  -I*o23
        G23     0		-G13	0		0		0		0	   -I*o23	I*o23
        I*o12	-I*o12	 0       -L12	0		I*o23	0		0		0
       -I*o12	I*o12	 0       0		-L12'	0		-I*o23	0		0
        0       0		 0		I*o23	0		-L13	0		-I*o12	0
        0       0		 0		0		-I*o23	0		-L13'	0		I*o12
        0       I*o23	-I*o23	0		0		-I*o12	0		-L23	0
        0      -I*o23	 I*o23	0		0		0		I*o12	0		-L23'
    ];

    n=expm(M*t)*n0;
end